Numerical study of the stability regions for half-quantum vortices 

in superconducting Sr2Ru04 
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Abstract 

We numerically solve the coupled Landau- Ginzburg-Maxwell equations for a model of a spin 
triplet p x + ipy superconductor in which whole or half-quanta of flux thread through a hole. We 
recover the pattern of stable and unstable regions for the half-flux quanta observed in a recent 
experiment. 

PACS numbers: 74.20.Dc,74.20.Rp,74.25.Ha,74.70.Pq 



I. INTRODUCTION 

Superconductors with triplet p x + ip y pairing are interesting because they can host half- 
quantum vortices with Majorana core states and non-Abelian braid statistics [TH3]. There 
is indirect evidence that the layered superconductor Sr2Ru04 has this pairing [I], and a 
search is on for "smoking gun" signatures that will confirm this. One signature -- chiral 
edge currents — has proved elusive, but recent work j5] has found striking results suggesting 
that half-quantum vortices have been detected. If this result is correct, it strongly supports 
the spin triplet pairing character of the superconducting order parameter. 

In the experiment reported in jS] a micron-sized annular flake of Sr 2 Ru04 is mounted on 
a cantilever. Its magnetic moment is monitored as magnetic fields both perpendicular (B x ) 
and parallel {B z ) to the c-axis are applied. As B z is increased, moment jumps corresponding 
to the entry of single-flux- quantum vortices into the hole in the annulus are easily observed 
(see figure [I]) . 
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FIG. 1: SRO sample and magnetization curves showing integer flux transitions adapted from [5]. 



When a sufficiently large in-plane field B x is applied, these entry-event jumps break up into 
two separate events, each with one-half of the original magnetic-moment jump. (See Figure 

B) 

The most obvious explanation is that we are seeing half-quantum vortices: The large 
B field presumably has broken a spin orbit coupling that had held the spin-triplet order- 
parameter d vector parallel to the pair angular momentum vector 1. The d vector can now 
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FIG. 2: Magnetization curves for various values of applied in-plane field. Adapted from [5]. The 
shaded region highlights the wedge shaped character of the half-flux state's stability region. 

rotate freely in the x-z plane and this freedom permits the existence of a half-quantum vortex 
defect around which the d vector and the order-parameter phase <fi both rotate through an 
angle ir while leaving the spin-triplet order parameter single valued. These rotations have 
the effect of producing a phase-winding vortex in one spin component, while leaving the 
other with no phase winding. A circulating "anti-vortex" current is nonetheless induced in 
the second spin component by the magnetic field, and at large distance the spin-up and spin- 
down flow velocities become equal and opposite. This means that there is a long-range spin 
current surrounding the vortex. It is the necessity of reducing the logarithmically divergent 
free-energy-cost of this spin current that mandates the use of small annular samples [6] . 

As B x is increased, the separation between the half-quantum jumps becomes larger. 
Vakaryuk and Leggett [7j have proposed that this phenomenon can be explained by a kine- 



mafic spin polarization. The different flow velocities for the up and down spin condensates 
give rise to an analogue of the Bernoulli effect which increases the magnitude of the order 
parameter for one condensate and decreases it for the other. The resulting magnetic polar- 
ization lies parallel to B x , and so alters the energy cost for a vortex to enter the sample. 
Assessing whether or not the energy gain from Vakaryuk-Leggett mechanism is sufficient 
to explain the experimental data requires a detailed accounting of the energy in the three 
dimensional magnetic field surrounding the sample and in the degree of polarization of the 
condensate. In this paper we perform this accounting by obtaining numerical solutions to an 
appropriate set of coupled Maxwell-Landau-Ginzburg-equations. With the specific geometry 
of the samples used in the experiment, and with reasonable values of the Landau-Ginzburg 
parameters, we find that we are able to qualitively reproduce the experimental data. 

In section |TT] we assemble the ingredients of the two-component Landau-Ginzburg free- 



energy functional that we use to model the superconductivity. In section III we explain how 
we find the field configurations that minimize our free energy functional, and how we extract 
from them magnetization curves that we can compare with the experimental data. Section 



IV presents the numerical results, and, finally, section M contains a brief discussion of these 



results and their significance. 

II. TWO-COMPONENT LANDAU-GINZBURG FORMALISM 

A. Spin and charge stiffness 

For a spin-triplet p x + ip y superconductor with a fixed orbital angular momentum vector 
1 = z, the order parameter is matrix- valued and of the form 

-(di + id 2 ) d 3 
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(1) 



The d vector has unit length, and for our application we will assume that it is perpendicular 
to a spin-quantization axis e 3 , which need not be the z axis. We therefore set d% = 0, 
d\ + id 2 = e^, and define the phases of A^ and A^ to be 9^ + ir and 0± respectively. Then 

X = \(0, + 9 l ), (2) 

0=^t~^)- (3) 



The authors of [6] write the free-energy density in the London form 
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This expression contains only the Goldstone fields x an d 4> and so ignores the free-energy 
cost of gradients in the magnitude of the order parameter. Nonetheless Q captures the 
essential far-from-core vortex energetics. In particular, in a half-quantum vortex either 9^ 
or 9_i (but not both) rotate through ±27r. Then x an d 4> rotate through ±7r. Far from the 
vortex core the B field will adjust so as to make the first term in K zero, so the total flux 
threading the half-vortex is given by 
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The remaining term cannot be screened by the B field and gives a contribution to the vortex 
energy that is logarithmically divergent at large distance. This divergent energy cost means 
that a half-quantum vortex be disfavoured unless the spin stiffness p sp i n is small and a finite 
size to the superconducting region cuts-off the long-distance contribution. 

The angle-valued Goldstone fields are not suitable for numerical work as they are not 
singled- valued in the presence of vortices. We need to write the free energy in terms single- 
valued fields. We therefore introduce fields ifo = \ip^\ exp{iQ^} and ip± = \ipi\ exp{z#j,}. The 
simplest form for the Landau- Ginzburg free-energy density that has the correct symmetries 
contains the kinetic-energy density 



K } 



h 2 



landau 



2m* 



V 



2ie 



A Vt 



V 



2ie 



A ^ 



26J t ■ J 



where 



ieh 
m* 
ieh 
m* 



r t h- 2 -^A]^-^(v + 2 -^A]r t 



n 



Pi v 



2ie 
2ie 



2ie 



(6) 

(7) 
(8) 



The current-current interaction introduces no new magnitude-gradient free-energy cost, and 
so does not affect the coherence length. 

If we set ip^ = \if)\e f, etc., and temporarily ignore derivatives of the common magnitude 
l^l , then 
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which is to be compared with desired London form of [6J. As p sp i n must be positive, b must 
be less than unity. Being just less than unity encourages half-quantum vortices. 

The presence the term 2bJ^ ■ Jj, in the free energy density alters the current J maX weii that 
couples to the magnetic field. We have 

8e 2 b 
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B. Kinematic spin polarization 

Legget and Vakaryuk propose [7j that the dependence of the half quantum vortex stability 
region on the applied in-plane field B x can be understood via the existence of a spontaneous 
spin polarization in the half quantum vortex state that arises from the difference between 
the spin-up and spin-down condensate velocities v-)- and v^ in the half-quantum vortex. 

The polarization occurs because for each of the two spin components in the Landau- 
Ginzburg theory we have a version of Bernoulli's equation: 

-m |v u | + = const. (11) 
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Here p u = \ip^,i\, and 



u(p) = ap + -f3p 2 , (12) 



is the potential part of the Landau- Ginzburg free energy density. 

As a consequence, the faster the superflow, the lower the order-parameter density. The 
polarizing tendency is proportional to 

l v tl 2 -lnl 2 = ( v t + n) • ( v t - n) ( 13 ) 

" * charge ' *spin v / 

This polarizing tendency gives rise to spin magnetic moment 

Here ps denotes the Bohr magneton, and g is a phenomenological parameter that we expect 
to be of order unity. (If IV'Ifi. were the actual density of cooper pairs, and if each of the two 



electrons in the pair contributes a Dirac moment of fieiectmn ~ /•*#, we would have g = 2.) 
The induced moment couples to the magnetic field to give a free-energy contribution 

AF = -B ■ ^ spin . (16) 

Jang et al. [5] assume that this moment lies in the x-y plane and so it is affected only by 
B x . We therefore account for it by including a term 

AF = -g f i B (\^\ 2 -\^\ 2 )\B ll \ (17) 

in our free-energy functional. 

C. Anisotropy 

The superconductor Sr 2 Ru04 is quite anisotropic, and it is necessary to replace the usual 
Landau-Ginzburg scalar mass m* with a mass tensor M* so that the free energy becomes 

Fy>,r,A] = Jj 3 x{j ( v + x A ) r ' (Mrl ' ( v ~ ir A ) ^ 

+«h/f + f H 4 + ^|V x A - B cxt | 2 | . (18) 

We choose a coordinate system such that the crystal ah planes lie in the xy-plane and are 
displaced from each other in the ^-direction. Then, the effective mass tensor takes the form 

M* = diag(m||, m||,mj_) (19) 

where the ratio 
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can be calculated from the measured penetration depths A^ and Ay of Sr 2 Ru04 and is found 
to be approximately 400. 

III. NUMERICAL METHOD 

Our goal is to create equilibrium magnetization curves that can be directly compared 
with those in [5] (see figures IT] and pq) . Because of the asymmetric and three dimensional 
character of the experimental samples, these curves have to be found numerically. 



For the numerical calculation, it is convenient to express the Landau-Ginzburg free energy 
in dimensionless form as 



F[Vt,^,A]= [ d 3 x\y2 
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Here J i = Re{— £0?(« -1 V — iA)tp*} is the dimensionless current, and b and \x are dimen- 
sionless parameters corresponding to the current-current coupling & in (J6J) and g/is in (17), 
respectively. The parameter k ~ 2.3 is the ratio of the in-plane penetration depth to the 
in-plane coherence length. Also 



M; ed = diag(l,l,ml/m;). 



(22) 



To find the magnetic field and condensate configurations that minimize the free energy 



functional (21) we used the commercial finite-element method solver COMSOL. 
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FIG. 3: On the left is shown a representation of the simulated geometry. On the right is a top-down 
view of the simulated ring. The axes in both figures are labelled in units of An. 



We will focus on the results for the single annular sample shown in figure [3j The ring's 
inner radius is 2 A|| while the outer radius averages 5.29 A|| where A||(T = 0) is given as 
152 nm for S^RuO^ in [I]. The height of the sample is 1.7 A||. This geometry is designed 



to approximate the experimental sample shown in Figure 111. The ring is centered in a 
cylindrical volume Q of height 20 Am and radius 20 A||. The dimensions of the cylinder were 
chosen to balance the competing effects of increased computation time for larger cylinders 
against spuriously high magnetic field energies caused by confining the field in too small a 
volume. 

The externally imposed magnetic field B = (B x , B y ,B z ) is established by imposing the 
inhomogeneous Dirichlet boundary condition 

A x = ^{.B y z-B z y), 
A y = ~(B z x- B x z), 
A z = -(B x y-B y x). (23) 

on the surface the dQ of the cylinder. 



The boundary conditions we impose on the order-parameter fields in (21) are the "natu- 
ral" boundary conditions that arise from the variational problem of free energy minimization. 
In other words we require the vanishing of the integrated out variation terms on the surface 
of the superconductor. This leads to 

n .(M: cd )- 1 .(--tA-^J i )^ = 

n ■ (M^)- 1 • ( ^ - iA - §J t ) fa = (24) 

on the surface of the superconducting ring. 

We find the fields that minimize the free energy by first choosing initial conditions for 
the magnetic vector potential and for the order-parameter fields, and then allowing them 
to relax to equilibrium. As the applied magnetic field along the z-axis is increased, we 
see vortices enter the superconductor and increase winding numbers about the ring. Each 
entry results in a discontinuous step in the magnetization. These transitions turn out to be 
difficult to model reliably. Although vortices enter the sample, the regions of metastabilty 
are large and geometry dependent. We therefore adopted a strategy that is similar to the 
field-cooling used in the actual experiments. We begin by imposing initial conditions 

fa = exp[m t 0]; fa = exp[in±<f>] (25) 



10 



that correspond to a selected flux state. (Here, is the azimuthal angle around the ring. For 
instance, (n^, %) = (0, 0) corresponds to the zero flux state while (1, 0) and (0, 1) correspond 
to half flux states.) After selecting the desired winding numbers, and setting the externally 
imposed magnetic field, we allow the system to relax to a local minimum of the free energy. 



This free energy is then calculated by evaluating the integral in (21). In this way we are 



able to construct diagrams of free energy versus applied magnetic field in the z-direction 
for different flux states and for different values of in-plane magnetic field. These diagrams 
reveal which flux state is energetically favourable at each value of applied z-axis field. Using 
these plots, we numerically compute the derivatives of the free energy with respect to the 
applied z-axis field, and thus construct the magnetization curves. 

IV. RESULTS 
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FIG. 4: This figure shows an example of calculated relative free energies of the integer and half- flux 
states near the first transition. Each red line represents the free energy of the half-flux state at 
different values of applied in-plane magnetic field. Higher magnetic fields are seen to result in 
relatively lower free energies of the half-flux state. 



Figure 0] illustrates the free energy of the system versus applied z-axis field (from to 
32.5 Gauss), and for a few values of applied in-plane magnetic field. The blue line depicts 
the (0,0) and (1,1) integer-flux states while the red line depicts the (1,0) half-flux state. [9] 
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One can see from the diagram that the half-flux state has a reduced free energy versus the 
integer-flux states for higher values of in-plane field. 
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FIG. 5: Computed magnetization curves obtained using p sp /p s = 0.25 and a magnetic dipole 
moment of 0.55//£. 



Applying a parabolic fit to the free energy curves, one can take the the derivative with 
respect to the applied magnetic field to obtain magnetization curves. Figure [5] displays 
magnetization curves versus z-axis field for various values of in-plane field as in Figure 2. 
The magnetization curves are seen to be qualitatively similar. By adjusting b and jl one may 
obtain the desired minimum in-plane stabilization field and stability region growth rate. 

The dimensionless current-current coupling parameter b can be related to the ratio of the 

spin and charge fluid densities via 

ps P _ 1 -b 
p s 1 + b ' 



(26) 
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Typical values of b which matched the experimental data were 0.4-0.6 corresponding to 
Psp/ps between 0.25 and 0.43. 

If p is interpreted as the magnetic dipole moment per particle of the spin condensate, it 
gives a magnetic dipole moment in units of \/2e 2 B c \ 2 m~ 1 . Using B c = 230 Oe and A = 152 
nm, the value of this unit is approximately 2.12 x 10~ 23 J/T w 2.2/xg. Typical values of the 
magnetic dipole moment which fit the data were found to be around 0.6 ps- 

The magnitude of the in-plane spin magnetic moment for the case of 240 Oe in-plane field 
is plotted in Figure [61 Near zero applied z-axis field the system is in the (n^, n±) = (0, 0) 
state. The presence of the in-plane magnetic field induces an in-plane magnetic moment 
even in the absence of kinematic spin polarization. Thus, there is a non-zero moment even in 
the integer flux states. The presence of kinematic spin polarization accounts for the sudden 
increase in moment at the transition to the half-flux state near 12 Oe. This additional 
moment vanishes when the system transitions to the unit fluxoid state at 16 Oe. 

V. DISCUSSION 

We have shown that the model described in section II can qualitatively reproduce the 
results of the experiment described in [5] . Using kinematic spin polarization and a current 
coupling term in a two-component Ginzburg-Landau model, numerical results show that an 
in-plane field results in the increasing stability of a half-quantum vortex state. This leads 
to a wedge-like stability region in the half-flux state versus applied c-axis field as shown in 
Figure [5j Moreover, this is achieved using physically reasonable values of p and p sp /p s - 

In comparing Figure [5] to Figure [2] there are some disrepencies that should be acknowl- 
edged. Firstly, the periodicity of the numerical result is approximately 30 Oe while that of 
the data is 16 Oe. This is due to the fact that, assuming a penetration depth of 152 nm, the 
simulated ring's hole diameter is 0.6 pm. while the actual sample has a hole diameter of 0.75 
pm. A smaller hole necessitates a larger applied field in order to achieve the flux necessary 
to induce a fluxoid transition. 

A more perplexing difference is in the magnitude of the moments, the numerical result's 
moments being an order of magnitude larger than those shown in Figure [2] These small 
magentic moments seem to be found only in the smaller samples examined by Jang, et. al. 
The moment data collected from larger samples, as in Figure 2 of [8], seems to be of the 
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FIG. 6: The calculated in-plane spin magnetic moment for 200 Oe of applied in-plane field. 

appropriate order of magnitude. At the present time, we can only speculate that the small 
moments are due to weakening of the superconductivity caused by crystal defects incurred 
during fabrication of the rings. 

This discrepency in the magnitudes of the magnetic moments makes the interpretation 
of the experiments in [5] difficult. One possible scenario explaining the half-height jumps in 
magnetization were Abrikosov vortices piercing the side-wall nearly half-way between the top 
and bottom of the ring. One of the several arguments against this was that the magnitude 
of the induced spin magnetic moment \ini is an order of magnitude less than that produced 
by a side- wall vortex. The estimates of (Xhi were determined by Jang, et. al. by using the 
formula [ihi = 5H z Afi z /4:(H x — H XtTn i n ) where A/z 2 is the jump in magnetic moment upon 
entry of a unit vortex, SH Z is the width of the stability wedge, and H x — H x ^ min is the amount 



14 



of applied in-plane field over the minimum necessary to see half-flux states. For the data 
shown in Figure [2j the implied hhi ~ 9 x 10~ 15 emu while, from Figure J6j \xm ~ 3.5 x 10 -14 
emu. So, since the measured A/j, z and hhi appear to be poorly understood this argument 
loses some credibility. A better understanding will require an understanding of the character 
and stability of wall-vortex states. 
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